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We introduce a model of attractive penetrable spheres by adding a short range attractive square 
00 ' well outside a penetrable core, and we provide a detailed analysis of structural and thermodynamical 

properties in one dimension using the exact impenetrable counterpart as a starting point. The 
■ model is expected to describe star polymers in regimes of good and moderate solvent under dilute 

conditions. We derive the exact coefficients of a low density expansion up to second order for the 
S_i , radial distribution function and up to fourth order in the virial expansion. These exact results are 

0- used as a benchmark to test the reliability of approximate theories (Percus-Yevick and hypernetted 

chain). Notwithstanding the lack of an exact solution for arbitrary densities, our results are expected 
to be rather precise within a wide range of temperatures and densities. A detailed analysis of some 
limiting cases is carried out. In particular we provide a complete solution of the sticky penetrable- 
sphere model in one dimension up to the same order in density. The issue of Ruelle's thermodynamics 
stability is analyzed and the region of a well defined thermodynamic limit is identified. 
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S ' I. INTRODUCTION 

Unlike simple liquids, where two-body potentials describe interactions at the atomistic levels, complex liquid inter- 
actions are always a result of an average process over the microscopic degrees of freedom. As a result, highly simplified 
models often accurately describe a number of experimental features ranging from structural to thermophysical prop- 
erties. Examples include colloidal dispersions, macromolecules, and combinations of the two A noteworthy feature 
of these system is that the hard-core repulsive barrier for very short range is not an essential ingredient of the model. 
In the case of highly ramified polymers in good solvent (star polymers), for instance, the center-of-mass of two polymer 
chains can be at a distance much smaller than their respective radii of gyration and they are well described by an 
effective Gaussian interaction ■ The simplest of this class of minimal bounded potential is the so-called penetrable 
£f~) • spheres (PS) [|[ that has attracted considerable attention in the last few years (see e.g. Ref. [|J and references therein). 
In this case the infinite barrier of the hard sphere (HS) potential is replaced by a finite one, thus allowing for a finite 
probability of penetrating inside the core. 
55 ' A major advantage of the PS potential is, of course, simplicity. On the other hand, it lacks an attractive part 
which is expected to be relevant in such a complex environment in view of the ubiquity of van der Waals dispersion 
forces. The purpose of the present work is to address this point by proposing a variation of the PS model in which a 
square well (SW) is added outside the core. This model, hereafter referred to as penetrable square- well (PSW), has 
an extremely rich phenomenology notwithstanding its simplicity, including a number of interesting limiting cases as 
will be discussed later on. 

One-dimensional bounded interactions do not belong to the class of nearest-neighbor fluids for which the total 
potential energy can be written as 

N-l 

U N {xi, . . . ,X N ) = ^ _ x i\) > (!•!) 



where cj>(r) is the pair potential and {xi, i = 1, . . . , N} are the coordinates of the N particles confined in a segment of 
length L, which eventually may be let to go to infinity. A necessary (but not sufficient) condition for a one-dimensional 
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fluid to satisfy Eq. (jl.ip is to be a hard-core fluid, i.e., a fluid made of particles which cannot penetrate one another 
due to the existence of an infinite repulsive potential barrier in <j)(r). 

Nearest-neighbor fluids admit an analytic exact statistical-mechanical solution in one-dimension [f|: the partition 
function, equation of state, and correlation functions of any order can be calculated analytically from the knowledge 
of the pair potential. In addition to HS [H, 0, 0, H S Gi|, both structural and thermophysical properties can be 
analytically obtained in one dimension for Baxter's sticky hard-sphere (SHS) potential [ill, fl^ . and for the SW 
potential [131 ] , but the technique permits in principle the analysis of a large class of nearest- neighbor one-dimensional 
potentials. 

In the absence of the nearest-neighborhood constraint (as happens with bounded potentials), the situation is far 
more complex and we are not aware of any general analytical approach to the problem, even in one dimension. As 
a matter of fact there exist only a few examples of analytically solvable one-dimensional models of this type, which 
include the Kac potential [lj] and the Coulomb potential [HI, For PS, it was observed [1, [I?) that the exact 
analytic solution for HS can be efficiently exploited to build a rather precise, albeit approximate, solution of the 
penetrable counterpart. This analysis is here extended to PSW interactions. Using a low density expansion and the 
corresponding exact solution for the SW problem, we derive the exact result up to the second order in a density 
expansion of the radial distribution function and up to fourth order in the virial expansion of the equation of state. 
These exact low density calculations arc contrasted with approximate theories such as the Percus-Yevick (PY) and 
the hypernetted chain (HNC) closures, thus providing an assessment on the relative reliability of both approximations 
and the low density expansion. As a preliminary simplified step to our calculation, we also examine the penetrable 
counterpart of the SHS problem, denoted as sticky penetrable spheres (SPS) in the following, which provides a 
guideline to tackle the more difficult PSW problem. 

The introduction of an attractive part of the potential into a penetrable interaction raises the important issue of 
the existence of a well defined thermodynamic limit [H,[l{|. We address this problem for the PSW model and provide 
compelling arguments to identify the stability region, which is guaranteed for a sufficiently small (« 0.5) ratio between 
the attractive and repulsive energy scales and arbitrary values of the other parameters. 

The remaining of the paper is structured as follows. In Section|TT]we introduce the model along with all its limiting 
cases (including the SPS fluid) and we study its stability. Section ITTT1 briefly accounts for the main equations necessary 
for the analytical solution of the nearest-neighbor class of fluids with arbitrary interactions. The exact solution of the 
one-dimensional SHS potential is derived within this general approach in Section IIIIBI and this is used to obtain the 
corresponding low density solution of the SPS in Section ITVl A similar analysis is carried out in Section [V] for PSW 
and the results are contrasted with those stemming from PY and HNC closures. Section IVIII contains some closing 
remarks, whereas some of the more technical details are confined in suitable appendices. 



II. THE PENETRABLE SQUARE- WELL (PSW) MODEL 

The penetrable square-well (PSW) model is defined by the following pair potential (see Fig. [T]) 

!+e r , r < a, 
-e a , cr<r<a + A 7 (2.1) 
0, r>a + A 7 

where e r and e a are two positive constants accounting for the repulsive and attractive parts of the potential, respec- 
tively. Here a is the diameter of the sphere (length of the rod in one dimension) and A < a is the width of the 
well. This model has a number of relevant limiting cases. When e r — » oo it reduces to a square-well (SW) fluid, 
whereas e a — > yields the penetrable-sphere (PS) model studied in Ref. Q in the one-dimensional case. In addition 
it gives rise to an interesting variation, referred to as "sticky penetrable spheres" (SPS), within an appropriate limit 
of infinite well depth of vanishing width (see below). Finally, we recover the hard-sphere (HS) fluid in the combined 
limit e r — > oo and e a — > 0. 

It is worthwile noting that the PSW model (and its variants) considered here is different from other app arently 
similar models like the Widom-Rowlinson's of interpenetrating spheres [20], the concentric-shell model \2l\ , or the 



pcrmcablc-sphcrc model [22 1 



As usual, a very important role is played by the Mayer function 

f(r) = e-*»-l, (2.2) 



where f3 = l/fc^T is the inverse of the thermal energy (ks is the Boltzmann constant and T is the absolute tempera- 
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FIG. 1: (Color online) Sketch of the penetrable square-well (PSW) potential (left column). The right column shows a few 
limiting cases: the sticky penetrable spheres (SPS) potential (e Q — > oo and A — * 0), the penetrable square-shoulder (PSS) 
potential (e r > —e a > 0), and the hollow hard-spheres (HPS) potential (e r = and e a < 0). 



ture). In the present model, this becomes 

!-7r, r <cr, 
+7a, o < r < a + A, 
0, r>cr + A, 

= 7r/ HS (r) + 7 a[0(r - a) - 0(r - a - A)} , (2.3) 

where 

7, = 1 - e"^ (2.4) 

is the parameter measuring the degree of penetrability varying between (free penetrability) and 1 (impenetrability) 
and 

la = e' 3 ^ - 1 (2.5) 

plays a similar role for the attractive part. Here fns( r ) = — cr) — 1 is the Mayer function for a HS model, which 
can then be recovered in the limit 7 r — > 1 and j a — > 0, and 0(r) is the usual step function equal to 1 for r > and 
otherwise. It also proves convenient to introduce the ratio 7 = 7 a /7r which is a measure of the depth of the attractive 
well, relative to the "penetrability" of the core. In that way, Eq. (|2.3[) can be rewritten as 

f(r) = 7r {/hs(t) + l[0(r — a) — 8(r - cr - A)]} 

= 7r/swM, (2.6) 

where fsw( r ) is the Mayer function of a SW fluid with the change 7 a — > 7. 

Although the PSW model can be defined for any dimensionality of the system, throughout the remainder of this 
paper we will specialize to the one-dimensional case. 
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TABLE I: Summary of the models. 



Model 



Acronym 



Penetrable Spheres 
Penetrable Square- Well 
Sticky Penetrable Spheres 
Penetrable Square-Shoulder 
Hollow Penetrable Spheres 
Hollow Hard Spheres 



PS 
PSW 
SPS 
PSS 
HPS 
HHS 
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A. The thermodynamic stability issue 

As anticipated in the Introduction, in our model we need to make sure that the system is always stable in the 
sense that the total energy is always bounded from below by —NB, N being the number of particles and B being an 
arbitrary positive constant p^ |. The physical origin of such an instability can be traced back to the fact that a soft 
core allows the possibility of a "collapsed-state" where the energy is no longer proportional to the number of particles 
N and a well-defined thermodynamic limit may not exist. In a classic paper, Fisher and Ruelle [l8l | provided a set of 
conditions on the pair potentials which are sufficient for stability, but the actual implementation of such conditions 
in soft-core systems is far from being trivial as was recently shown for Gaussian-core models (23j and Lennard- Jones 
fluids [H. 

In the PSW model, the issue is clearly related to the interplay between the two energy scales, e r and e a , for the 
repulsive and attractive parts of the potential. As shown in Appendix [A| we predict that the system might be unstable 
when e r < 2e a whereas we prove that it is certainly stable in the opposite case e r > 2e a . 



B. The sticky limit: the SPS model 



It is instructive at this point to consider a particular limit of the PSW model which will be referred to as the sticky 
penetrable-sphere (SPS) model (see Fig. [1]). This is a variation of a widely used sticky hard-sphere (SETS) model 
introduced long time ago by Baxter [25[,which has proven to be extremely useful in the framework of complex fluids, 
recently even in its anisotropic version [261 ] . The simplest way of introducing it is at the level of the Mayer function 
[see Eq. ([2.3p ] , which becomes 



where 



/shs (r) = /hs (r) + aad + (r - a) 

e(r)-9(r-a) 

6+{r) = hm , 

a^o+ a 



The relation with the SPS model is then provided by 



/sps (r) = 7r/sHS (r) 



(2.7) 



(2.8) 



(2.9) 



In the original SHS model [25|], a = 1/12t > 0, (r playing the role of an effective temperature) but the connection 
with the PSW model is readily achieved from Eq. (j2.6[) by considering the limits A — > and e r — > oo so that 
a = (7 /7 r )(A/cr) remains finite. In spite of its usefulness, the SHS model is known to suffer from some mathematical 
drawbacks, the most important of them being that it is unstable in spatial dimensions greater than 1, as pointed out 
by Stell [27| , in view of the divergence of the virial coefficient corresponding to a close packed configuration. For the 
SPS model we will be able to achieve a number of exact results which can be exploited as a guideline for the more 
complex PSW model. 



C. Other limiting cases 



In all previous cases, we have tacitly assumed e a > 0. In principle, however, nothing prevents to consider the 
opposite case e a < (what implies j a < 0). In this case the PSW potential gives rise to an interesting class of 
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models, at least from an academic point of view, with two positive energy scales (e r and |e a |). If e r > e a |, we get a 
purely repulsive potential that could be called the penetrable square-shoulder (PSS) model (see Fig. [TJ. A peculiar 
situation occurs if e r < \e a \: when two particles approach they have first to overcome the barrier \e a \ at r = a + A: 
once done, they feel an attractive well of depth |e a | — e r for r < a. Thus the potential is attractive for short distances 
and repulsive for larger distances. The simplest version of models with e a < and e r < \e a \ corresponds to e r = 0, 
which will be referred to as the hollow penetrable-sphere (HPS) model (see Fig. [T]). If, in addition, the limit |e | — > oo 
is taken, one gets an athermal potential that will be referred to as the hollow hard-sphere (HHS) model since the 
particles look like hard spheres of diameter a + A with a "hole" of diameter a inside. If two particles arc separated 
by a distance larger than a + A, they behave as hard spheres and never know about their holes. On the other hand, if 
the separation between them is smaller than a, they can never separate a distance larger than a. In the HHS model, 
7 a — > — 1 and j r — > 0, so that the functions 2/2 ( r ) and g 2 (r) are well defined (see below). In Section [VI] we will discuss 
the results for representative values of the parameters. A summary of the penetrable models treated in this paper, 
along with the corresponding values for the values for e r , e a , and A characterizing them, is reported in Table [J 



III. BASIC FORMALISM FOR EXACT PROPERTIES OF NEAREST-NEIGHBOR POTENTIALS 



A. General scheme 



The great advantage of dealing with one-dimensional models is that they are usually amenable to exact solutions, 
at least in the limit of sufficiently short range interactions [l4| ■ The trade-off is of course the fact that these models 
do not have phase transitions. In the context of fluids, this translates into the fact that there exist exact solutions for 
the HS, SHS, and SW models [f| 0, H, HH, EH, HH, H§] . The same formalism allows one to tackle non-nearest- neighbor 
one-dimensional fluids [i| thus leading to an approximate solution. Let us recall the main results of this approach, 
referring to Ref. [l3T | for a self-contained treatment. The main quantity to be computed is the Laplace transform of 
the Boltzmann factor e _/3 ^ < - r - ) : 



fi(s) 



ds 



-W(r) 



This is directly related to the Laplace transform of the radial distribution function <?(r), 

/■OO 

G (s) = / ds e~ sr g (r) , 



where p = N/L is the density of the one-dimensional fluid. The relation is (see Refs. [10|>[l3| for details) 

n -1 



G(s) 



1 



Here £ is a solution of the equation 



dlnfi(s) 



- f 



(3.1) 



(3.2) 



(3.3) 



(3.4) 



Finally, the equation of state (EOS) (and hence the whole thermodynamics) can be cast into the very simple form 



(3.5) 



where P is the pressure. 

In practice, the scheme goes as follows. Evaluate fi(s) from the Boltzmann factor by a Laplace transform, Eq. (|3.1|) ; 
solve for £ from Eq. (|3.4[k insert the result into Eq. (|3.3p : invert the Laplace transform (|3.2p to obtain g(r) and, in 
parallel, compute the EOS from (|3.5p . 

As a final remark, we anticipate that, when dealing with discontinuous potentials (or Boltzmann factors), it is 
convenient to introduce the cavity function y(r) which is related to the radial distribution function g(r) and the pair 
potential <p(r) by the general relation 



g (r) = e-nWy{r) 



(3.6) 
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Moreover it can be expanded in powers of the density 

oo 

y(r) = 1 + J2p n yn(r). (3.7) 

>i=i 

In principle, the knowledge of all y„ coefficients provides the exact solution to the cavity function y(r) (provided 
that the above series converges) and hence to the problem. This also allows us to assess the reliability of well known 
approximations involving the direct correlation function c(r) and the cavity function [3(3| , such as the Percus-Yevick 
(PY) closure 

c(r) = f(r)y(r) (3.8) 

and the hypernctted chain (HNC) closure 

c(r) = f(r)y(r) + y(r) - 1 - lny(r). (3.9) 



B. Exact solution of the SHS model in one dimension 



Let us particularize the above procedure to derive the exact solution of Baxter's SHS model in one dimension. 
Starting from the Boltzmann factor 



its Laplace transform (|3.ip yields 



e -P4>(r) = Q(r- a )+ao-6 + (r-a), 



O(s) = I aa + - I e 



Equation (|3 .4[) can then be arranged to get the following quadratic equation 

£, 2 a 2 a (1 - pa) + £er (1 - pa) - pa = 0. 

Its physical solution is 

^1 + 4apa/(l - pa) - 1 



which can be substituted into Eq. (|3 . 5[) to give 



2aa 



(3.10) 



(3.11) 



(3.12) 



(3.13) 



(3.14) 



PP _ v/l + 4apcr/(l - pa) - 1 
p 2apa 

which represents the EOS for this system. In order to get the exact radial distribution function, we exploit Eq. (|3.3 
to get 



G 



1 OO 



E 



(3.15) 



We can now use the binomial theorem to expand (aa + -^j) n and invert the Laplace transform (|3 . 2|) term by term 
by using the residue theorem, to obtain 



g( r ) = ^ ip n (r ~ na)Q(r - na), 



n=l 

1 1 



p(aa+ l/i) r 



(aaT5 + (r)+J2[l)(av) 



k=l 



i-k 



r k—l e —£r 



(3.16) 
(3.17) 



which is the correct result found in Ref. [1 lj | with a different method. 
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IV. EXACT PROPERTIES OF THE SPS MODEL 

Next we turn our attention to the corresponding penetrable SPS counterpart. Following Ref. the basic idea 
hinges on deducing the exact low-density orders of the SPS model from those of the SHS, which can be evaluated 
exactly. Each term y n (r) can be represented as a sum of suitable diagrams, whose forms for y\(r) and yi (r) were 
given in Ref. [|| and will not be repeated here. Each bond in the diagrams corresponds to a Mayer function /(r) and 
the ones for the SPS and SHS models are related by Eq. (|2.9[) . For the SHS model previously discussed, the exact 
cavity function y(r) does not have a Dirac delta function at r = a and the regular part is continuous at that point. 
None of these two properties is any longer true for the SPS model, as further elaborated below. Here and in the 
following we set a = 1 for simplicity. 

The result is 



Vi 



(SPS) 



r) 



2 (SHS) , \ 



(SPS) , n 3 (SHS) , s , 4 (SHS) M . 7 r (SHS) / , . 7r (SHS) , n 

v\ [r) = %v\a ( r ) + 2 7 r y 2 B ( r ) + yy 2 c ( r ) + y^n ( r ) 



(4.1) 
(4.2) 



where the first-order density term is 



y (SHS) (r) = ( 2 _ r _2a)9(2 



r) + a 2 [25+ (r) + 5+ (r - 2)] . 



(4.3) 



Note that this has a delta singularity at r 

[- (3 - r 2 



(SHS) / n 



and is continuous at r 
+ 6a(l-a)] 6(1 - r) 



1. For the second order in density we have 



(3 



-a 3 [3<5 4 



3a (3- 
f 5+ (r 



[6(3-r) -0(l-r)] 



3)] 



(4.4) 



(SHS) / x 
V2B ( r ) 



— (6 — 2r — r 2 ) — a (6 — 6a — r) 



9(l-r) 



-(2-r) (4-r) 



4a - 3r) 



[0(2-r)-e(l-r)] 



(SHS) / \ 

v\c ( r ) 

(SHS)/ \ 
V2D W 



+a 2 (1 - 2a) [25+ (r) + 5+ (r - 2)] - 2a 3 5+ (r - 1) . 
[^ SHS) (r)] 2 , 

[-(3 - 2r) + 2a(3 - 3a - r)]9(l - r) 

+ [-(2 - r) 2 + 4a(2 - a - r)][9(2 - r) - 9(1 - r)} 

+2a 3 [65+(r) + 5+{r - 1) + 35+(r - 2)] - a 4 [4(5 2 (r) 



2)] • 



The functions (|4.4[) — (|4.7[) present some peculiar properties. In particular, (i) the regular parts of J/ 2 A HS ^( r )' V2B 
and 2/2i) iS ' > ( r ) are discontinuous at r = 1; (ii) y^ 1 ^^)! 2/2S HS ''( r )' an< ^ y^zj 13 ^ 7 ") nave a delta singularity at r = 1 
and (iii) y^?^ i r ) ari< i 2/2 S D iS '' ( r ) present delta-square singularities at r 



(4.5) 
(4.6) 



(4.7) 

(SHS) / 



'(r), 



and r = 2. However, these three classes 



of singularities cancel out when setting j r = 1 in Eq. 



to obtain the total second-order function y 2 SHS ''( r ) [3l| - 

On the other hand, since for SPS y^t* ( r )> 2/2S HS ''( r )' 2/2C iS ^( r )' an< ^ V2D HS \ r ) are weighed by different powers of 
7 r (3, 4, 4, and 5, respectively), the corresponding exact second-order cavity function is discontinuous at r = 1 and 
has a delta singularity at r — 1 and delta-square singularities at r — and r = 2. The delta singularity at r — 1 is 
responsible for a diverging fourth virial coefficient of the SPS model (see below). 



V. EXACT PROPERTIES OF THE PSW MODEL 



A. Calculation of yi and j/2 



As already mentioned, the SPS model suffers from the same drawbacks as the original SHS model plus some 
additional ones, so that it can hardly be regarded as a sound model in higher dimensions. However it has served 
as a test bench for analytical techniques. Armed by these tools we can now tackle the more difficult PSW model 
which has the SW fluid as a reference model. We recall that the latter does not have an exact solution in higher 
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dimensions but it is amenable to an exact treatment in one dimension [13]. The discussion follows closely the route 
already introduced for the SPS model, namely the density expansion, Eq. (|3.7p . The radial distribution function g(r) 
is related to the cavity function y(r) by Eq. ()3.6|) which with the help of Eqs. (|2.2p and (|2.3|) yields 



(! - 7r) V (r) , r < 1, 
g (r) = { (1 + 77r ) y (r) , 1< r < 1 + A, 
y(r), r > 1 + A. 



(5.1) 



As in the SPS model, the cavity function can be exactly computed up to second order in density, this time by reducing 
the problem to the solution of the SW model. 
The first order term reads (A < 1) 



yi (r) = 7r < 



2(1 + 7 2 A) -r (I + 27 + 27 2 ) , < r < A, 

2 - 27A - r, A < r < 2, 

7 (2 + 7) (r - 2) - 27A, 2 < r < 2 + A, 

(2 + 2A-r) 7 2 , 2 + A<r<2 + 2A, 

0, 2 + 2A < r. 



(5.2) 



The second order can be reduced to the calculation of the corresponding diagrams of the SW model as anticipated. 
We find 



V2 (r) = 7^ W) (r) + 2 7 ^T (r) + |,f c W) (r) + f ^ (r) 



(5.3) 



where the explicit calculation of the various terms is described in Appendix [Bl and is given by Eqs. (|B1[) . (|B2[) . (|B4|) . 
and (|B8|> . It can be checked that these expressions reduce to those of the SPS model, Eqs. (|4.4p - (|4.7|) . in the limit 
7 — » 00 and A — > with a = 7A/ a = const. 



B. Computation of B2, Bz, and B 4 

The EOS can be obtained from the knowledge of the radial distribution function g(r) through a number of routes. 
The most common ones arc the virial route 



BP f°° r) 

f— = Z(p,p) = l + 2 d ~\ d p j dr r d y (r) — / (r) 



(5.4) 



the compressibility route 



and the energy route 



dp 



X (p,(3) = l + 2 d dv d p dr r*- 1 [g (r) - I] 



(5.5) 



U - I m d 
- = u{p,(3) = - 



l + 2 d v d pf3 / dr r d - x 4>{r)g(r) 



(5.6) 



where d is the dimensionality of the system and v d = {ir/4) d / 2 /T(l + d/2) is the volume of a d-dimensional sphere of 
unit diameter. Thermodynamic consistency for the exact g(r) requires the three routes to be completely equivalent 
and hence 



x _1 (p,/3) = i-\pz{pM 

dp 



(5.7) 
(5.8) 



For an approximate g(r), on the other hand, the consistency is no longer granted and different routes (or combinations 
of them) may lead to different results. 
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Let us specialize to the one-dimensional case of the PSW model, where we have just derived the exact g(r) up to 
second order in a density expansion. Equations (|5.4[) - (|5.6[) become, using the potential (|2.1[) . 



Z(p,/3) = l + P7 r[(l+7)y(l)-7(l + A)y(l + A)], 



X {p,P) = l + 2p{ / dr [(l- 7r )y(r)-l] + 



l+A 



dr [(1 + ~1rl)y{r) - 1] + 



+ OC 



dr [y (r) - 1] 



l + A 



e r (1 - j r ) / dr y (r) - e a (1 + 7r7) 



l+A 



dr y (r) 



Inserting the expansion (|3.7p for the cavity function y(r) we find 

Z = 1 + B 2 p + B 3 p 2 + B iP 3 

X = 



1 + X2P + X3P 2 + XiP 3 + 



1 

2(3 



+ U 2 p + U 3 p 2 + M 4 p 3 + 



(5.9) 
(5.10) 

(5.11) 



(5.12) 
(5.13) 

(5.14) 



Clearly Eq. (|5.12p is the virial expansion for the compressibility factor Z whereas (|5.13[) and (|5.14p are the analogous 
expansions for the isothermal compressibility x an d the energy per particle u. If the exact coefficients y n appearing 
in Eq. (|3.7p are known, the above three quantities provide the identical exact EOS. 

On starting from the second-order values B 2 , X2, and u 2 one can obtain perturbativcly higher orders term by term 
from the knowledge of yn{r)- The result can be cast into the form 

(5.15) 
(5.16) 

n>3, (5.17) 
(5.18) 

Note that B n depends upon y n -2 so that knowledge of the exact y\ and y 2 allows the computation of the exact virial 
coefficients up to B4. The third- and fourth-order results can be obtained from Eqs. (|5.2p and (|5.3p . After some 
algebra, one gets 

B 3 = 7 3 [1 - 7 A (2 - A - 2 7 A)] , 
Xa = 4B 2 2 -3B 3 , 

^7 r (l + 7,7)A(2-A-4 7 A) 



B 2 


= Ir 


(1- 


7A) , X2 = 


-2B 2 , 


u 2 = e r (1 - 7r) - e a (1 + 7r7) A, 


B n 


= Ir 


[(1- 


1-7) Vn-2 (1) - 


-7(14 


- A) y n -2 (1 + A)] , n>3, 


Xn 


= 2 


(1- 


-7) / dry n - 


2(r) + 


(l + 7r7) / dry n - 2 (r)+ / dry n ^ 2 (r) 

Jl Jl+A 


u„ 


= (r 


(1- 


7 r) / dry n . 
Jo 


2 (r) 


/■l+A 

e a (1 + 7r7) J dr 2/„_ 2 (r) , n > 3. 



(5.19) 
(5.20) 



"3 = 


e r 2 ft 

(1 - 

1 9 R 


7r ) [3 - 2 7 A (2 - 




2 dp 3 ' 




Bi = 


7 6 r 

2 L 


- 7 A (3 - 3A - 6 




7 5 r 
+ lH 7 - 


- 7 A (21 - 15A - 




2 L 


- 7 A(12-6A- 


X4 = 


-4 (2B 3 2 
1 9 R 


— 3B 2 B 3 + B4) , 


U4 = 


3^^ 4 - 





(5.21) 



.,3 A2^ 



- A 2 + 3 7 A 2 + 3 7 2 A 2 - 3 7 3 A 2 )] , (5.22) 

(5.23) 
(5.24) 

The three routes provide consistently identical results for B 3 and B4, i.e., the relations (|5.7p and (|5.8p are verified, as 
they should. In the energy case the following identity is needed: 

d <9 7r d #7 d 

d . e a (1 + 7 r 7) - e r 7 (1 - 7r ) 5 



£r (1 - 7r) « H 



7r 



97 



(5.25) 
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Equation (|5 . 22[) gives the exact fourth virial coefficient as a function of the three relevant parameters of the PSW 
model, namely j r , 7, and A. The results for the PS and SW models are recovered as 

}™ B * = (TW = * (4 + IT - 2 ) ' (5 ' 26) 



lim B A = 1 - 7A ( 3 - 3A - 67A + -A 2 + -7A 2 + 5 7 2 A 2 , (5.27) 

e r ^oo y 2 2 J 

respectively. On the other hand, while B 2 and B3 arc well defined in the SPS limit (7 — > 00 and A — ► with 
a = 7A = finite) [see Eqs. (|5 . 1 5[) and (|5.19p ]. the presence of the terms 7 3 A 2 in Eq. (|5.22p implies that B4 — > 00 
in the SPS model. Equation (|5.16[) shows that this is a direct consequence of the divergence of y^ PS \r) at r = 1. 

However, 2/2 (1) = finite, so that B4 is well defined in the SHS model (jy r = 1), as shown by Eq. (|5.27|) . 

The second, third, and fourth virial coefficients for the PSS model (e a < 0) arc still given by Eqs. (|5.15[) . (|5.19[) . 
and (|5.22p . except that 7 < 0. In the case of the HPS model (e a < and e r — > or, cquivalently, 7 a < 0, 7 = ja/'Jr, 
and 7 r — > 0), one gets 

lim B 2 = - 7a A, lim B 3 = 0, lim B 4 = -j^A 3 . (5.28) 
The special case of the HHS model is obtained by further taking the limit e a — > 00 ("/ a — > — 1). 



VI. SOME ILLUSTRATIVE CASES AND COMPARISON WITH THE PY AND HNC 

APPROXIMATIONS 

The approximate character of a given closure can be typically inferred by looking at g 2 (f) along with the corre- 
sponding fourth virial coefficient -B4. Being coefficients in a density expansion, both can be either positive or negative. 
We now plot the exact g 2 {r) and B\ for some illustrative cases and compare them with the PY and HNC theories 
(JHIBD and (|33]) . The PY and HNC results corresponding to Eq. (jO]) are [3(j 

VV (r) - iHT (r) + HvfP (r) , (6.1) 



y 2 HNC (r) = 7^ W) (r) + 2^.yg W) (r) + \yf^ (r) . (6.2) 

Comparison with Eq. (|5.3p shows that the HNC theory neglects y^n ( r ) and the PY theory neglects, in addition, 

2/2C? ( r )- ^ s a consequence, the expression for the fourth virial coefficient in the PY and HNC approximations 
depends on the thermodynamic route. The corresponding results can be found in Appendix [Q 

Let us start with g 2 (r). As a prototypical PSW system we have chosen e a /e r = 0.25 and A = 0.5. Figure [2] 
shows gz(r) for ^T/e r = 0.5 and fc^T '/e r = 1. It can be observed that the HNC and PY tend to overestimate and 
underestimate, respectively, the values of (72 (f) in the overlapping region r < 1. This is due to the fact that y^D^ ( r ) 
is generally negative in the region r < 1, while y^^\r) is positive definite and larger than the magnitude of y^^\r). 
Inside the well (1 < r < 1 + A) the PY and HNC curves practically coincide at fcsT/e r = 0.5, both being rather 
inaccurate, while at the higher temperature ksT/e,. = 1 the PY prediction is quite good. Moreover, the PY theory 

is a better approximation than the HNC theory for r > 2. This is explained by the fact that y^c^ ( r ) ^V-n^^i 7 ") = 
in the region r > 2, so that g\ Y (r) coincides with the exact g 2 (r) for r > 2 in the case of the SW model (7,. = 1). If 
7 r < 1 the combination y^c^M + IrV^D ( r ) docs not vanish for r > 2 but is still rather small for the cases of Fig. 

[21 For r > 2 + 2A, both 5 2 1NC ( 7 ') and ff2 Y ( r ) become exact since y^c^'W an d y^D / \ r ) vanish in that region. 

Figure [3] depicts the function 32 (^) for a representative case of the PSS model (see Section Til C[) . Most of the 
preceding comments in connection with Fig. [2] apply here as well. Finally, the function g 2 (r) corresponding to the 
HPS model is shown in Fig. 2] for A = 0.5 and two temperatures: fcsT/|e a | = and fcsT/|e a | = 0.5. Note that the 
zero-temperature case is equivalent to the HHS limit. It is interesting to note that the curves corresponding to both 
temperatures are quite similar, except for a change of scale. In the HPS model the HNC theory gives the exact gz(r) 
because, for large I7I, y^^ir) scales with 7 4 , while it has a weight 7^ and so does not contributes to y 2 (r). Similarly, 
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FIG. 2: (Color online) Plot of the second-order radial distribution function, <?2(r), for a PSW model with e a /e r = 0.25, A = 0.5, 
and fcsT/e r = 0.5 (top panel) and kBT/e r = 1 (bottom panel). The solid, long-dashed, and short-dashed lines correspond to 
the exact result, the HNC approximation, and the PY approximation, respectively. 



y-n^\ r ) scales with 7 3 and so it does not contribute to y^ir) either. On the other hand, JrV^^^i 7 ") 1S different from 
zero in the regions < r < A and 2 < r < 2 + 2A and it is there where the PY theory fails, yielding Y {r) = 0. 

In order to have a feeling of the behavior of the exact B4, we now plot them for some representative values of the 
parameters. Figure [S] shows the exact [see Eq. (|5.22p ] and the approximate (see Appendix (Cj values of the fourth 
virial coefficient as functions of temperature for the same PSW model as considered in Fig. [21 i.e., the one defined by 
fa/fr = 0.25 and A = 0.5. While the exact B4 goes to —00 as T — > 0, the HNC and PY theories artificially predict 
a divergence to +00. We can observe that the best agreement with the exact curve corresponds to B 4 ' c up to 
ksT /e r ~ 0.5 and to B^ Y,c thereafter. The worst behaviors correspond to Bf ,v = i?^ NC,e and B^ Y ' V . 

The results for the PSS model considered in Fig. [3l namely e a /e r = —0.5 and A = 0.5, are displayed in Fig. [6] 
For low temperatures this model reduces to the HS model of diameter 1 + A. It is found that B V A Y ' v and, especially, 

PY c 

B A ' present an excellent agreement with the exact B4. On the other hand, the poorest performances arc presented 

, dHNC.d D HNC.e , D PY,e 

by B A ' = B A and B 4 

We have also evaluated B4 for the HPS model at various values of fcsT/|e a |, as depicted in Fig. [71 and compared 
with the PY approximation (compressibility route). As said before, the HNC theory becomes exact for the HPS 
model. Interestingly, in this case both the virial and the energy routes from the PY approximation yield exact results, 
even though y\ Y (r) is not exact. 

It is worthwhile noting that B4 is not a monotonic function of temperature in the PSW model (see Fig. [5]): it is 
negative for low temperatures, reaches a positive maximum value at an intermediate temperature, and then decays, 
reaching a very small negative minimum value at a certain temperature, and finally going to zero from below. Although 
hardly apparent in Fig. [SI the behavior of B4 is also non-monotonic in the PSS model: it is generally positive and 
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FIG. 3: (Color online) Plot of the second-order radial distribution function, g2(r), for a PSS model with e a /e r — —0.5, A = 0.5, 
and fcsT/e r = 0.5 (top panel) and fcsT/e r = 1 (bottom panel). The solid, long-dashed, and short-dashed lines correspond to 
the exact result, the HNC approximation, and the PY approximation, respectively. 



decays as the temperature increases, but eventually reaches a very small negative minimum value and thereafter tends 
to zero from below. In contrast, the fourth virial coefficient of the HPS model (see Fig. [7]) is negative definite and 
monotonically increases with increasing temperature. 



VII. CONCLUSIONS AND OUTLOOK 



In this paper, we have introduced the PSW model and outlined a number of exact results for this model in one 
dimension. The potential contains two energy scales (the core barrier e r and the well depth e a ) and two length scales 
(the core diameter a and the well width A). This model is a variation of the widely used square- well one with a finite 
energy barrier replacing the hard core. As such, this is not a nearest-neighbor system and there exists no general 
approach leading to an exact solution even in the one-dimensional case. In spite of this we have been able to obtain 
the exact first few coefficients in the density expansions of the relevant structural and thcrmodynamical properties. 
Specifically, we have computed both the cavity and radial distribution functions up to second order in density and the 
virial expansion up to fourth order. As a stringent test of the calculations, we have explicitly checked that different 
routes to thermodynamics (virial, compressibility, and energy) are consistent one another up to this order. 

This model includes a variety of other models as special cases. By taking the limit of infinitely narrow and deep 
well (e — > oo, A — > 0) we obtain the SPS model which can be also reckoned as a variant of the SHS model with 
penetrable core. Upon reversing the sign of the attractive energy scale e a we obtain a PSS model with successive 
soft repulsive barriers of decreasing height. If the second barrier is higher that the first (— e a > e r > 0), we find a 
potential that is attractive for short distances and repulsive for larger distances. An interesting situation, that we 
have denoted as HPS, corresponds to — e a > e r — 0. In the limit of zero temperature (or, equivalently, e a -co) 
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12 3 4 



r 

FIG. 4: (Color online) Plot of the second-order radial distribution function, gi(r), for the HPS model with A = 0.5 and 
kBT/\e a \ = 0.5 (bottom panel) and fesT/|e a | = (top panel, corresponding to the HHS model). The solid, long-dashed, and 
short-dashed lines correspond to the exact result, the HNC approximation, and the PY approximation, respectively. Note that 
the HNC approximation provides the exact result in the HPS model. 




FIG. 5: (Color online) Plot of the fourth virial coefficient, B4, as a function of fcjT '/e r for a PSW model with 6 a /e r = 0.25 
and A = 0.5. The solid, dashed, and dotted lines correspond to the exact result, the HNC approximation (virial/energy and 
compressibility routes), and the PY approximation (virial, compressibility, and energy routes), respectively. 
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Exact 



- -a- - HNC(v,e). 

HNC(c) . 

...a-py(v) 




0.0 0.5 1.0 1.5 

km 

B r 

FIG. 6: (Color online) Plot of the fourth virial coefficient, £>4, as a function of ksT/er for a PSS model with e a /er = —0.5 
and A = 0.5. The solid, dashed, and dotted lines correspond to the exact result, the HNC approximation (virial/energy and 
compressibility routes), and the PY approximation (virial, compressibility, and energy routes), respectively. 



0.00 




FIG. 7: (Color online) Plot of the fourth virial coefficient, B4, as a function of ksT '/\e a \ for a HPS model with A = 0.5. The 
dashed line corresponds to the exact result whereas the dotted line correspond to the PY approximation (compressibility route). 



the HPS model becomes the HHS model, characterized by an infinitely high barrier between a and a + A. Here the 
equilibrium state consists of "chains" of connected particles: two adjacent particles of the same chain move freely, 
provided that the distance between their centers does not exceed c; on the other hand, particles of different chains 
behave as hard spheres of diameter a + A. In the limit |e tt | — > (and also if — e a = e r > 0) the PSW fluid reduces 
to the PS one, and all results obtained here are consistent with previous analysis on the PS model within this limit. 
Finally, all results smoothly converge to the HS limit when e r — > 00 and \e a \ — ► (or — e a = e r — > 00), as expected. 

The combined effect of the absence of a hard core and the presence of a finite attractive part of the PSW potential 
raises the issue of the existence of a well-defined thermodynamic limit of the system. We have analyzed this issue in 
detail and we have assessed the limits of stability as a function of the ratio between the attractive and repulsive energy 
scales: when e a /e r < 5 the system is stable whereas in the opposite case the system might be unstable (when e r > 
and e Q < the system is always stable independently of the energy scales). The SPS limit turns out to be always 
unstable since the exact fourth-order virial coefficient diverges, unlike the corresponding SHS counterpart which is 
well-behaving. 
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A main advantage of exact relations is that one can asses the reliability of approximate theories. A comparison 
with PY and HNC closures unveils the corresponding strengths and weaknesses of both. We have found that each 
of them has a domain in space where it outperforms the other, and we have explained why this is so in terms of 
the exact and approximate behavior of the second-order cavity function 1/2 M- As a general feature, HNC tends to 
overestimate the cavity function within the core whereas PY has the opposite tendency. On the other hand, PY is 
consistently superior in the large-r region. Both approximate theories produce artifacts in the low temperature region 
of the fourth virial coefficients. 

It would be extremely interesting to extend the present work in some respects. While our analysis has provided a 
careful comparison of the PY and HNC approximate theories with respect to the exact result, we have not attempted 
a detailed physical interpretation of the results. This is because our exact analysis was limited to the lowest orders in 
density, which are expected to be valid only within a rather limited region of the phase diagram. It turns out, however, 
that even this limited knowledge can be exploited to construct rather precise approximations for the PSW model in 
the limits of small penetrability (1 — j r <C 1) and high penetrability (7,- <C 1), for arbitrary value of the density. This 
analysis mirrors that already performed for the PS model 0, [X^] , can be tested against numerical simulations, and 
nicely complements the exact low-density results presented here. The boundedness of the class of penetrable-sphere 
potentials raises the possibility of a phase transition even in a one-dimensional system (32l . [33| and the presence of 
the attractive part might also give rise to additional transitions in the fluid phase. We plan to address this point in 
future work. 

It is worth stressing that the large number of parameters present in the PSW fluid (two energy scales, two char- 
acteristic lengths, density, and temperature) may render the phase diagram analysis quite problematic, so an exact 
understanding of the low-dimensional behavior, where the analysis can be carried out almost fully analytically, is 
always welcome. Having done this, the extension to three dimensions should be facilitated, and our results predict 
an extremely rich phase diagram which might be useful to describe complex fluids with soft cores within a unified 
picture. 

APPENDIX A: RUELLE'S STABILITY CRITERION 

According to Ruelle's criterion, a sufficient condition of thermodynamic stability is [IH, [l!| 

iV-l N 

U N (x u ...,x N )=J2 4>{\xi-Xj\)>-NB (Al) 

i=l j=i+l 

for all configurations {xi}, where B is a fixed bound. 

Let us first demonstrate that for small repulsion (e r < 2e a ) there exists at least one configuration violating the 
stability constraint. We consider a particular configuration where the N particles are distributed into N/s overlapping 
clusters of s particles each, so that the centers of the s particles belonging to the same cluster coincide (with a tolerance 
A/2). Next, the centers of two adjacent clusters are separated by a distance between a and a + A. The potential 
energy corresponding to this configuration is 

„ , s Ns(s-T) [N \ a , 

Un{s) = L -=— ^ - - - 1 s 2 e a , A2 

s 2 \ s J 

The first term on the right-hand side represents the repulsive energy of the s(s — l)/2 pairs of each cluster, times the 
number of clusters. The second term is the attractive energy of the interaction between the s particles of each cluster 
and the s particles of its neighbor cluster, times the number of pairs of adjacent clusters. Obviously, the value of the 
total potential energy Un depends on the value of s. The extreme cases are s = 1 and s = N. We then see that the 
value that minimizes Un(s) is 

which is meaningful only if e r < 2e a . The corresponding minimum value of Un(s) is 

2 



U N (s,) = -N 



€ - + ^e a 2-^ 
2 16 V e a 



e r < 2e a . (A4) 



The quantity enclosed between brackets grows linearly with N and so it is not bounded. Therefore, if e r < 2e a , 
there exists at least a configuration that violates Ruelle's criterion. On the other hand, we note that if e r > 2e a the 



16 



minimum of Un(s) is reached at s = 1, in which case 

U N (s) > U N (s = 1) = — (JV - l)e„ , e r > 2e a , (A5) 

so that all these special configurations are consistent with Ruelle's criterion. Indeed we now show that no other 
configurations violate Ruelle's criterion if e r > 2e a so that the model is thermodynamically stable if the above 
condition is satisfied. 

Without loss of generality we can see any given configuration of N particles as a set of M clusters (1 < M < N), 
each cluster i being made of Si overlapping particles (i.e., any pair of particles of a given cluster are separated a 
distance smaller than a). For fixed M and {si}, the total potential energy can be decomposed as 

U N ({ Si }; M) = t#*"({*}; M) + U% tei ({si}; M), (A6) 

where 

M 

2 



U^({ Si };M) = |^ Si ( Si -l) (A7) 



! = 1 



is the contribution associated with pairs of particles inside each cluster and U^ tcl is the contribution associated with 
pairs of particles belonging in different clusters. Note that in the latter contribution the energy for each pair can be 
e r (if the separation is smaller than a), — e a (if the separation lies between a and a + A), or zero (if the separation is 
larger than a + A). It is clear that the minimum value of U N lteI is achieved when all the particles of a cluster interact 
attractively with all the particles of the neighbor cluster: 

M-l M-l 

UN tei ({si};M) > -e a s » s i+i > - £ « s l s l+1 -e aSl s M . (A8) 



Therefore, 



M M 



U N ({ Sl }; M) > s ^ -!)-yE s '( Si -! + Si +i) = WN({si}\ M) , (A9) 

j=l i=l 

where so = sm and sm+i = si- Given M, what is the set of population numbers {si} that minimizes Wn subject to 
the constraint that J2i=i s i = Using a Lagrange multiplier A, the problem reduces to solve 



_d_ 

dSj 



M 



W N {{ Si };M)-\J2si 



= ^(2s j -l)-e a (s j _ 1 + s j+1 )-X = , l<j<M. (A10) 



The solution is s, = N/M and A = (N/M)(e r — 2e a ) — e r /2. This could have been expected by symmetry arguments. 
Therefore, given M clusters, the minimum Wn is obtained with a uniform distribution Sj = s = N/M. For e r > 2e a 
we can thus write 

U N ({ Si }; M) > W N ({ Si = N/M}; M) = ^M(N/M)(N/M - 1) - e a M(N/M) 2 

= (e r /2 - e a )N 2 /M - Ne r /2 > -Ne r /2 , (All) 
which proves that Ruelle's stability criterion is satisfied. 



APPENDIX B: SECOND-ORDER CAVITY FUNCTIONS FOR SW 

The first-order term y^ W \r) (for A < a = 1) is given by Eq. (|5.2[) with j r = 1. This allows for a straightforward 
determination of y^^^ i r ) as 



y2c W) ( r ) = yf' ( r ) 



(SW) 



(Bl) 
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Next, one can also evaluate the Fourier transform of the integral corresponding to the 2 A diagram. Going back to 
real space, the result is 

viTir) = -^7(l + 7) 2 (l-A-r) 2 e(l-A-r) + ^(l+ 7 )(l + 2 7 + 3 7 2 )(l-r) 2 e(l-r) 

(2 + 4 7 + 3 7 2 ) (1 + A - rf 6 (1 + A - r) + ^ 7 2 (1 + 7 ) (1 + 2A - rf 6 (1 + 2A - r) 
-\ (1 + 7) 3 (3 - rf 9 (3 - r) + ^ 7 (1 + 7 ) 2 (3 + A - rf (3 + A - r) 

(1 + 7) (3 + 2A - rf 6 (3 + 2A - r) + ^ 7 3 (3 + 3A - rf 9 (3 + 3A - r) . (B2) 

For y^B ( r ) wc can ma ke use of the identity 

/oo 
ds yi (s)f(s)f(\r-s\), (B3) 
-oo 

which leads to the result 

2/E W) W = 7 (1 + 7) (2 - A - 4 7 A) (A - r) 9 (A - r) + 7 (1 + 7 ) 2 (1 - A - rf 9 (1 - A - r) 

- (1 + 3 7 + 5 7 2 + 3 7 3 ) (1 - rf 9 (1 - r) + 7 (2 + 4 7 + 3 7 2 ) (1 + A - rf 9 (1 + A - r) 

-7 2 (1 + 7 ) (1 + 2A - rf 9 (1 + 2A - r) + i (1 + 7 ) 2 (4 - r - 4 7 A) (2 - r) 9 (2 - r) 

-7 (1 + 7 ) (4 - r - 4 7 A) (2 + A - r) 9 (2 + A - r) + ^ 7 2 (4 - r - 4 7 A) (2 + 2A - r) 9 (2 + 2A - r) . 

(B4) 

Computation of y^^ (r) is much more laborious and requires a different route. We go back to the general formalism 
and compute the exact f2(s) from the Laplace transform (|3.ip . which is 

fi(s) = lI(l + 7 - 7e - sA ). (B5) 
s 

Equation (|3.4|) then yields for the parameter £ the following density expansion: 

£ = p + (1 - 7 A) p 2 + [1 - 7 A (2 - A - 2 7 A)] p 3 + ■ • ■ . (B6) 

Inserting this solution into Eq. (|3 ,3[) and inverting the Laplace transform (|3.2[) , we can obtain the corresponding radial 
distribution function g 2 i r ) correct up to second order in density. Use of Eq. (|5.1|) then yields the corresponding cavity 
function (r) and then y^^\r) is given by the difference 

(r) = 2y2 (r) _ 2 ,f/) (r) _ 4y f/) (r) - 4 y f c w ) (r) . (B7) 

This provides the result for r > 1. Inside the core we have three different regions under the assumption that A < 1, 
namely 0<r<A,A<r<l — A and 1 — A < r < 1. The quadratic expression in each region can be obtained by 
imposing continuity conditions and with some help from numerical evaluation. The final analytic result is 

y^ w) ( r ) = -27 (1 + 7) [(1 + 7 + 7 2 ) r - 2 + A (l + 3 7 - 7 2 )] (A - r) 9 (A - r) 

-7 (1 + 7 ) 2 (1 - A - rf 9 (1 - A - r) + (l + 3 7 + 5 7 2 + 3 7 3 ) (1 - rf 9 (1 - r) 

-7 (2 + 4 7 + 3 7 2 ) (1 + A - rf 9 (1 + A - r) + 7 2 (1 + 7) (1 + 2A - rf 9 (1 + 2A - r) 

+ (1 + 7 ) 2 [r (1 - 2 7 - 7 2 ) - 2 + 4 7 + 2 7 2 + 4 7 A] (2 - r) 9 (2 - r) 

+4 7 2 (1 + 7) (r - 2 - A + 7A) (2 + A - r) 9 (2 + A - r) - 7 4 (2 + 2A - rf 9 (2 + 2A - r) . (B8) 

Note that the first derivative y'(r) is discontinuous at r = A, 2, 2 + A, and 2 + 2 A, as it can be inferred from its 
explicit computation at this order in density. 
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APPENDIX C: CALCULATION OF B 4 FOR THE PSW MODEL IN THE PY AND HNC 

APPROXIMATIONS 

Here the fourth virial coefficient predicted by the PY and HNC approximations from the various thermodynamic 
routes (v = virial, c = compressibility, e = energy) are given. 

1. PY approximation 

Using Eq. (|6.ip . along with the recursion relations (|5.15p ~ (|5.18[) . we have 

B FYv = f r [3 - 7 A (9 - 7A - I67A + A 2 + 67A 2 + 6 7 2 A 2 - 2 7 3 A 2 )] 

-Jr. [4 - 7 A (12 - 6A - 18 7 A + A 2 + 3 7 A 2 + 3 7 2 A 2 - 3 7 3 A 2 )] , (CI) 

X 4 PY = -4 (2Bl 3B 2 B, + B P 4 Y - C ) , (C2) 

B PY,c = Jjr ^ _ ^ _ _ + g A 2 + 167A 2 + lg7 2 A 2 _ 4 7 3 A 2~)j 

o 

-Jr. [4 _ 7 a (12 - 6A - 18 7 A + A 2 + 3 7 A 2 + 3 7 2 A 2 - 3 7 3 A 2 )] , (C3) 

O 

ul Y = lr(\ + 7r7 ) A [12 - 18 7r - 6A (1 - 2 7r + 6 7 - 10 7r7 ) + A 2 (l - 2 7r + 6 7 - 26 7r7 + 9 7 2 - 36 7r7 2 

-7 3 

-12 7 3 + 16 7r7 3 )l e a - Jr. (1 _ 7r ) A [16 - 28 7r - 6 7 A (6 - ll 7r ) + 6 7 A 2 (3 - 8 7r + 6 7 - 14 7r7 ) 
6 

-7 (1 + 7) A 3 (3 - 10 7r + 3 7 - 28 7r7 )] e r . (C4) 
The fourth virial coefficient associated with the energy route, B^ Y ' e , is obtained from Eq. (|C4[) as 

B 4 PY ' e =3fd/?'< Y (/J'). (C5) 
Jo 

PY e 

Its expression is quite long and so it is omitted here. In addition to its dependence on j r and 7, B 4 ' depends on 
e a /(n a e a - n r e r ) with (n„, n r ) = (1, 1), (1, 2), (1,3), (1, 4), (2, 1), (2, 3), (3, 1), (3, 2), (4, 1). 

It is instructive to consider some special cases. First, the results for the PS model correspond to the limit e a — > 
(7 -> 0) or e a -> -e r (7 -> -1): 

dPY,» 

lim5 4 PY ^ I™ -^-^=7^7,-2), (C6) 

^'-J^^-^^- 2 )- (C8) 

In the special case of the HPS model (e a < and e r — > 0) one finds that B^ Y ' V and -B4 Y ' e reduce to the exact result 
[see Eq. KM ] but lim er ^ B% Y '° = ~7 4 A 3 . 
The conventional SW model corresponds to e r — > 00 (7,. — > 1): 

lim B 4 PY '" = 1 - 7 A ^3 - 4A - 7 7 A + iA 2 + ^ 7 A 2 + ^ 7 2 A 2 - i 7 3 A 2 ") , (C9) 



V 2 2 ' 2 ' 2 
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lim Bl Y > c = 1 - 7 A ( 3 - 3A - 6 7 A + ^A 2 + ^ 7 A 2 + ^ 7 2 A 2 , (CIO) 

lim B 4 PY ' e = % - 7 A U - 3A - 6 7 A + ^A 2 + 5 7 A 2 + ^ 7 2 A 2 - J 7 3 A 2 ^ . (Cll) 
£ r ^oo 5 \ 2 2 2 J 

If, furthermore, the SHS limit ( 7 — > oo and A — > with 7 A = const) is taken in Eqs. (|C9[) - (|Clip . an artificial 
divergence of B4 is obtained. 

The results corresponding to HS are obtained by taking either the limit e r — > 00 ( 7r — > 1) in Eqs. (|C6[) - (|C8[) or the 
limit e a — > ( 7 — > 0) in Eqs. (|C9[) ~ (|C11|) . In either case one sees that the virial and compressibility routes yield the 
exact result, while the energy route value is wrong by a factor 4/5. A third possibility consists of taking the limit 
e a — > — 00 (7 — > —1) in Eqs. (|C9[) - (|C11[) . However, in this latter case the energy route yields an incorrect dependence 
on A: 

lim lim Bf Y ' e = % + — (6 + 6A + A 2 ) . (C12) 

The fact that the right-hand side of Eq. (|C12[) is not proportional to (1 + A) 3 implies that if one starts from Y ' e for 
the PSS model of shoulder height and width — e a and A, respectively, and then one takes the limit e a — ► — 00 to get 
the HS model of diameter 1 + A, the result has an artificial dependence on A. This anomaly of the PY description 
was discussed in Ref. HH . 

2. HNC approximation 

Similarly to the preceding analysis, from Eq. (|6.2|) one gets 

Bf Clt = |^T'°, (C13) 

= -4 (2Bl 3B 2 B 3 + B 4 HNC ^) , (C14) 

b hnc,c = ^ ^ 7 _ 7A (pi - 15A - 36 7 A + 3A 2 + 16 7 A 2 + 16 7 2 A 2 - 4 7 3 A 2 )] 

-JtL [4 - 7 A (12 - 6A - 18 7 A + A 2 + 3 7 A 2 + 3 7 2 A 2 - 3 7 3 A 2 )] , (C15) 

„HNC _ 1 <5 R HNC,« (n-\R\ 

u 4 -3QpB 4 , (C16) 

Equation (|C16[) implies that i? 4 4NC,e = B™ G ' V . This confirms that, in general, the energy and virial routes are 
thermodynamically consistent in the HNC approximation (35| . It is also noteworthy that the fourth virial coefficient 
predicted by the HNC approximation in the virial and energy routes is exactly three halves the one predicted by the 
PY approximation in the compressibility energy route, Eq. (|C13p . This simple relation is not restricted to ID models 
since it also occurs in the 3D PS model [l7|. It would be extremely interesting to check whether relation (|C13|) is a 
general property valid for any interaction potential and for any dimensionality. 
In the PS and SW limits Eq. (fCT5|) becomes 

respectively. 

The three routes in the HNC theory yield the exact result (|5.28[) in the HPS limit. However, as in the case of the 
PY theory, an artificial divergence of B4 is predicted in the SHS limit. 
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